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ABSTRACT 

We study the orbital structure of a self-consistent N-body equilibrium configuration of 
a barred galaxy constructed from cosmological initial conditions. The value of its spin 
parameter A is near the observed value of our Galaxy A = 0.22. We classify the orbits 
in regular and chaotic using a combination of two different methods and find 60% of 
them to be chaotic. We examine the phase space using projections of the 4D surfaces 
of section for test particles as well as for real N-body particles. The real particles are 
not uniformly distributed in the whole phase space but they avoid orbits that do not 
support the bar. We use frequency analysis for the regular orbits as well as for the 
chaotic ones to classify certain types of orbits of our self-consistent system. We find 
the main resonant orbits and their statistical weight in supporting the shape of the 
bar and we emphasize the role of weakly chaotic orbits in supporting the boxiness at 
the end of the bar. 

Key words: galaxiesistructure, kinematics and dynamics. 



1 INTRODUCTION 



The type of orbits that appear in various types of galaxies 
have been studied extensively up to now. (For a review 
see Contopoulos 2002). However relatively few studies 
have been made on the o rbital structure of N-bod y 
, systems (see for e xample ISparke Sellwoodl (Il987f) 
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ISundinet all (Il993ll. ISundinl (119961) . 
" 1998|1 . iBarnes and Tohlind (l200lh 
2OO4I ). ICeverino KlvpinI (|2005l )). In this paper 
consider a N-body model of a triaxial barred galaxy, 
spin parameter is relatively large, comparable to the spin 
parameter of our Galaxy. Our study aims to find the 
main types of orbits in such a galaxy. We emphasize the 
chaotic orbits that are 60% of the totality, in contrast with 
non-rotating N-body systems where the chaotic orbits are 
only 30%. The chaotic orbits occupy most of the space 
outside corotation. However the chaotic orbits are not 
completely mixed. They are separated for long times by 
cantori. Relatively few chaotic orbits escape from the 
system during a Hubble time. 

There are also many ordered orbits of various resonant 
types, mainly inside the bar. The topology of the orbits is 
governed by the resonances between the rotational and the 



* This work started by Dr. N. Voglis with Dr. M. Harsoula. It 
was very sad that Dr. Voghs passed away suddenly on 9/2/2007. 
After that Dr. G. Contopoulos continued this work and completed 
this paper in collaboration with Dr. M. Harsoula. 



epicyclic frequencies. The third dimension plays a minor role 
on the type of orbits. 

A comparison between the N-body orbits and the or- 
bits of test particles in the same potential shows that real 
orbits avoid systematically certain types of resonant orbits, 
occupying preferentially orbits that support the bar. In fact 
models consisting of particles uniformly distributed initially 
are not self-consistent. 

The paper is organized as follows: In section 2 we de- 
scribe the method used for constructing a realistic bar model 
and in section 3 we describe the characteristic features of the 
bar. Section 4 discusses the level of chaos. Section 5.1 de- 
scribes the various types of orbits by using surfaces of section 
and section 5.2 introduces a frequency analysis of the orbits. 
Finally in section 6 we summarize our conclusions. 



2 COSMOLOGICAL INITIAL CONDITIONS 
FOR CREATING A BAR-LIKE GALAXY 

For creating the final configuration of a bar-like galaxy we 
ha ve used two differ ent N-body codes: a tree code developed 
bv lHernauistldigSSll and a conservative technique code de- 
signed by Allen et aD (|l990h . Initially 5616 particles of equal 
mass are positioned in a cubic grid limited by a sphere and 
having a total mass equal to the mass of a galaxy.The system 
expands initially according to an Einstein-de Sitter model 
universe and then we impose a density perturbation with a 
spherical profile and a power law dependence on the radius 
r, i.e. 
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s(r) = ^ « I 

^ ' p ^(n+3)/2 



(1) 



Such a profile is consistent with a power-law spectrum 
of the density perturbation in the early Universe. Here we 
set n = — 2, a value that is in the range predicted by the 
Cold Dark Matter scenario for galactic mass scales (Davis et 
al. 1985). Furth er deta il s for these "quiet initial conditions" 
can be found in IVodisI (|l994l ). The unit length is Ikpc, and 
the unit of time tunit^lMyr while the total mass is taken as 
lO^^M©. Using the tree code, we find that the violent relax- 
ation process has been completed after t = 4000tunits. Then 
we split each particle into 27 particles arranged in a small cu- 
bic grid, while the sum of these children particles conserves 
the mass, energy and angular momentum of the parent par- 
ticle. The new configuration consists of about 1.4x10^ parti- 
cles and is now put as initial condition for the new run using 
the conservative technique code (hereafter c-t) for another 
150 half- mass crossing times. The c-t code gives a smooth 
potential containing 120 terms (mon opole, quadrupole and 
triaxia l) . All the details are given in Kalapotharako s et al.l 
In the same paper the time evolution of the potential 
coefficients is shown in figs. 5a,c. They seem to get stabilized 
fast and stay remarkably constant even for times compared 
to the Hubble time. The final output gives a triaxial ellipti- 
cal galaxy. Its ellipticity is about E7 on the Y-Z plane and 
about E5 on the X-Y plane, in a radius that corresponds to 
2 half mass radii (hereafter hmr). 

Since our aim is to create a system that is rotating 
around its smallest axis (i.e. the z-axis)we simply rotate the 
projections of the velocity vectors of all the particles on the 
X-Y plane to become perpendicular to their position vec- 
tors with a counterclockwise sense of rotation, and therefore 
we give the maximum possible rotatio n on this plane. Thi s 
mechanism has already been used in Voglis et al.l (|2006|). 
Then we let the system evolve using the c-t code again, for 
another 250 dynamical times and calculate the "spin param- 
eter" A (Peebles ,1969) as a function of the radius. The spin 
parameter is calculated by the following formula: 



J\E^ 



(2) 



where G is the gravitational constant. The angular momen- 
tum J, the bounding energy E (in the inertial frame), and 
the mass M are measured as cumulative quantities along 
cylinders on the X-Y plane, having the Z-axis as their main 
axis. 

At a time t = 250tunits the value of A has increased 
dramatically in the inner parts of the barred galaxy and 
stays close to A = 0.22 over most of the extent of the system. 
This value corresponds to the maximum value of A for a 
rotationally supported galaxy a nd is close to the one of ou r 
Galaxy, as it was calculated bv lEfstathiou JonesI (|l979h . 

In the c-t code we use a new runit that corresponds 
to the limiting radius of the whole bound system i.e. one 
runit is approximately 120 kpc in real units and a new tunit 
that corresponds to the half mass crossing time of the sys- 
tem (li^reafte^Junct^ time is ~ SOOhmct 
(see kalapotharakos et al.l (|2004[ ) for details), then Ihmct ^ 
45, GMyrs. 
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Figure 1. A snapshot of the isodensity contours at time=12/imct 
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Figure 2. The time evolution of flp of the bar, in radians per 
half mass crossing time(hmct) 



3 CHARACTERISTICS OF THE BAR 

In order to calculate the angular velocity of the pattern 
(hereafter Qp) of the system, we must find the time evo- 
lution of the angle between the great axis of inertia of the 
rotating system and a fixed axis, that is the Y-axis of the 
inertial frame of reference. However, by observing the time 
evolution of the isodensities on the rotating plane, we notice 
the existence of a differential rotation of the pattern, i.e. 
the inner parts seem to rotate faster that the outer parts 
at the beginning of the c-t run. This lasts up to a certain 
time, until a single Qp is established throughout the radius 
of the whole configuration. In Fig. 1 a snaphot of the iso- 
densities on the X-Y plane is shown at 12 half mass crossing 
times (hmct) from the beginning of the c-t run. We notice a 
trailing spiral structure, as the galaxy is rotating as a whole 
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counterclockwise. This spiral structure seems to travel out- 
wards before it vanishes later on. Therefore angular momen- 
tum is transferred to the material outwards. This is a well 
known mechanism that makes the bar grow stronger and 
slow down (Lynden-Bell and Kalnajs, 1972). These authors 
have shown that density waves in the form of temporary 
spiral structures carry angular momentum from the inner 
parts to the outer parts. This transference is due to the in- 
teraction of the bar perturbation of the gravitational field 
with stars moving in resonant orbits particularly near the 
ILR and corotation. They also proved that only in trailing 
spiral structures do the grav ity torques carry angular mo - 
men tum outwards. L ater on iTremaine &: Weinberg! (| 19841 ) 
and [Weinberg (|l985l l extended this analysis showing that 
angular momentum transference can be considered as a par- 
ticular type of dynamical friction. N-body simulations have 
confirmed this evolution scenario of bar-like galaxies (e.g. 
Athanassoula, 2003). 

In the present experiment this spiral structure vanishes 
after a short transient period of our calculations i.e. at t ~ 
20hmct. However, in experiments with a larger value of Qp 
these spiral structures can survive even for times comparable 
to a Hubble time (Voglis et al., 2006). 

In Fig. 2 we plot the time evolution of the angular ve- 
locity in units of rad/hmct of the material inside 2.0 rhm. 
The fluctuation of at the beginning of the c-t run is re- 
lated to the differential rotation of the pattern which trans- 
fers angular momentum outwards. This phenomenon creates 
spiral structures transiently (see Fig. 1) and causes the for- 
mation of tidal torques between the differentially rotating 
parts. Therefore the material inside the bar can temporar- 
ily accelerate or decelerate, before the whole system acquires 
the same Qp. The very small fluctuations after this transient 
time are less than 1^ /hmct. We notice that the material of 
the bar seems to slow down and stabilize its angular veloc- 
ity, after t ~ 150hmct to a value of Qp ~ 21" /hmct. This 
value of Qp is equivalent to {2ti) / {17 hmct)— 0.37 rad/hmct 
or to 8.25km/ {kpc. sec) (dashed line in Fig. 2). 

In order to study the orbital structure of this bar-like 
galaxy and the level of chaos we need good statistics and 
therefore a second multiplication of the particles is neces- 
sary. By considering only the bound part of our galaxy (in- 
side the radius 1.0 runit) we split again each particle into 
9 particles arranged in a small cubic grid around the parent 
particle, while the sum of these children particles conserves 
again the mass, energy and angular momentum of the par- 
ent particle. The new configuration consists now of about 
1.1x10^ particles. We then let the system relax for another 
50 half mass crossing times. After a short time period the 
system reaches again an equilibrium and rotates with the 
same Qp pattern: Qp = 0.37rad/hmct as a whole. 

For studying the orbits that support the bar and deter- 
mine the level of chaos we take a certain snapshot, where 
the bar is parallel along one of the main axes of the inertial 
frame of reference on the plane of rotation, freeze the po- 
tential of the system and calculate the orbits in the rotating 
frame of reference. 

At 50 tunits of the c-t run, after the new multiplication 
of the particles, the bar of the galaxy has been aligned with 
the Y-axis of the inertial system of reference. This is obvious 
in Fig. 3 where we plot the isodensity contours on the X-Y 
plane. 
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Figure 3. Isodensity contours on the X-Y plane at the end of the 
N-body run, together with the projections of the real particles of 
the galaxy. 

In Fig. 4 the surface density profile is plotted along 
the main axis of the bar (black line, Y-axis) and across the 
main axis of the bar (grey line, X-axis). It is obvious that an 
exponential- like profile can be fitted (thick line) in the first 
case almost until the limit of the bar at ~ 2hmr, which is a 
well known property, in the literature, for barred galaxies. 
The dashed line is the exponential fit of the surface density 
across the main axis of the bar. 

In Fig. 5 we plot the ellipticity=l — b/a, where a is the 
semi-major axis and b is the semi-minor axis as a function of 
the distance along the Y-axis. At this snapshot, the rotating 
frame of reference coincides with the inertial frame of ref- 
erence. The calculation has been done using the isodensity 
contours of Fig. 3. The ellipticity is close to E5 in the inner 
parts and drops to El near the end of the bar, while the 
isodensities become almost spherical outside corotation. 

The boxiness of the isophotes is a well known feature in 
barred galaxies. This effect can be measured by the shape 
parameter c, which can be determined from the equation of 
a generalized ellipse (Athanassoula et al. 1990): 

{\y\/ar + {\z\/br = l (3) 

where a and b are the semi- major and semi- minor axes (cal- 
culated from the equidensity contours) and c is the param- 
eter describing the shape of the generalized ellipse. For c=2 
we obtain a standard ellipse, for c > 2 the shape approaches 
a parallelogram and for c < 2 we have a lozenge. In Fig. 
6 we plot the shape parameter c as a function of the semi- 
major axis of each isodensity contour of Fig. 3. From this 
figure we conclude that in our model, the inner and the 
outer parts of the bar can be well fitted by ellipses while 
in a region in between we have somewhat orthogonal-like 
isophotes. This r esult agrees with the N-body experiments 
of I Athanassoula k MisiriotisI ()2002i ). 

An estimate of the perturbation of the bar is the value 
of Sq/qo, where is the mean density and is the devia- 
tion from it along an annulus in a certain distance from the 




4 A^. Voglis, M. Harsoula, and G. Contopoulos 




Y-AXIS / X-AXIS (half mass radii) 

Figure 4. The surface density profile at a time=250thmass along 
the Y-axis (black line) and along the X-axis (grey line). 
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Figure 6. The shape parameter c, as a function of Y (semimajor 
axis). 
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Figure 5. The ellipticity on the X-Y plane as a function of the 
semi-major axis of the bar. 
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Figure 7. The bar perturbation Sg/go as a function of y along 
the bar on the plane of rotation. The end of the bar is at ~ 2hmr. 



center. The calculations are made by projecting the parti- 
cles on the plane of rotation (x-y plane) and dividing the 
surface into annuli of equal width. Then we also split the 
xy plane in azimuth (with dO = 2 ^ 7r/40) into 40 successive 
angles creating bins. Then we count the particles in each bin 
to calculate the density p. We plot the maximum amplitude 
Sq/qo = {g — go) /go in every annulus as a function of y. By 
changing dO from dO = 2 * 7r/20 up to dO = 2 * 7r/50 we have 
found optimal values, such that a further decrease of dO does 
not change the maximum amplitudes. In Fig. 7 we present 
the maximum value of the perturbation of each annulus as 
a function of the radius. The maxima appear always along 
the bar (y-axis). We observe that the absolute maximum 
perturbation is close to 80% and is located a little inside 
the half mass radius of the system, while the value drops 
to ^ 25% near the end of the bar. For the region outside 



the bar we can only trust the mean value of the saw-tooth 
part. This perturbation is considered large and is often ob- 
served in real barr ed galaxies. See for example the paper of 
iLaurikainen et all (^004), where the m=2 and m=4 compo- 
nents are calculated for a sample of 112 spiral galaxies. 

Observers recognise the different components of a 
barred galaxy by plotting the density profile along and paral- 
lel to the the major or m inor axis of the bar as it is proposed 
bv lLutticke et al.1 (|2000l l. The bulge length is marked by the 
increased light distribution over the exponential disk well 
above the bar. In Fig. 8 we plot the surface density parallel 
to the minor axis for two different z. For z=0 (black line) we 
see that a central spherical component is present (we also 
mark the limits of the bar length on the x-axis). For z=1.6 
hmr a plateau appears that corresponds to the bar. There- 
fore we can conclude that this central outcrop corresponds 
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X-AXIS (half mass radii) 

Figure 8. Surface density along the bar minor axis for z=0 (black 
line) and z=1.6/imr (grey line). 
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Figure 9. The equipotential contours at the end of our calcula- 
tions. 



to a bulge that has been created self consistently when the 
system has reached an equilibrium. 

Another important quantity for the orbital study is the 
corotation radius. This can be determined from the maxi- 
mum value of the effective potential all along the main axis 
of the bar, as well as from the equipotential contours (Fig. 
9) where the Li and L2 points determine approximately the 
radius of corotation. Therefore we conclude that in our N- 
body galaxy, corotation is placed at approximately 2.2 hmr. 
The Langrangian points L4 and L5 are placed at approxi- 
mately 2.0 hmr. Figure 9 gives extremely smoothed contours 
as a result of the program used by " Mathematica" . 
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Figure 10. A snapshot of the orbits of the particles in the sample 
on the plane Log(AI) — Log(Lj) at 1200 radial periods. 



4 STUDY OF THE LEVEL OF CHAOS 

The purpose of this section is to classify the orbits of our 
self-consistent system in ordered and chaotic ones and find 
their energy distribution. 

This task can be accomplished by using a combination 
of two methods, that is the Specific Finite Time Lyapunov 
Characteristic Number (SFTLCN), or simply Lj, and the 
Smaller ALignment Index, (SALI), or simply Alignment In- 
dex (AI)(Skokos 2001, Voglis et al. 2002). 

The well known Finite Time Lyapunov Characteristic 
Number defined by the relation 

FTLCN =-log^-^^ (4) 

' \m\ 

is not very suitable for detecting chaotic orbits, because 
chaos is weak in the present case and the LCN of most of 
these orbits is small compared with the inverse of the dy- 
namical time of the system; therefore these chaotic orbits 
cannot be detected by integrating for times equal to only a 
Hubble time. 

This problem can be overcome by using the Specific 
Finite Time Lyapunov Characteristic Number (SFTLCN), 
or simply Lj for every orbit which is given by the formula: 



1=1 



(5) 



where tj is the integration time, Nj is the number of time 
steps (At = tj/Nj) and aij is the stretching number (Voglis 
and Contopoulos 1994) at the time step i {i = 1, ...Nj). The 
stretching number aij is defined by the equation 



In 
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\^j(U)\ 



(6) 



where ^j{ti) is the length of the deviation vector from the 



6 A^. Voglis, M. Harsoula, and G. Contopoulos 








Figure 11. (a) The energy distribution of the total number of particles, (b) The effective potential along the X,Y and Z axes, (c) The 
energy distribution of particles having chaotic orbits (black line) and of particles having regular orbits (grey line). 



orbit j at time t = ti. To avoid numerical overflows the 
deviation vectors are normalized to unity at regular time 
steps. The equations giving the deviation vectors are linear 
and therefore a change of their measure by a constant does 
not affect the information we get. The time evolution of the 
deviation vectors is found from the variational equations of 
motion. The coefficients of these equations are evaluated nu- 
merically. The quantity Trj is the average radial period of 
the orbit of the particle i.e. the time that a star needs to 
go from its pericenter to its apocenter and back. This period 
scales with the energy Ej of the orbit as Trj ~| Ej |~^/^ al- 
most independently of the angular momentum (Voglis 1994). 

The above definition of the SFTLCN allows us to ex- 
tend our calculations to different times for every orbit, up 
to a point when the integration time is long enough to allow 
LCNs to be stabilized at constant small values. 

In the second method (using AI) , we use the properties 
of the time evolution of the deviation vectors (Voglis et al. 
1998, 1999). 

In particular, we consider the time evolution of two 
arbitrary different initial deviation vectors and of 
the same orbit. If the orbit is chaotic, then the smallest 
of the two norms dj-{t) = \^ji{t) — ^j2{t)\ and dj^{t) = 
l^ji(t) + ^j2{t)\ is called Alignment Index (AI) and tends 
exponentially to zero. However, in our calculations we im- 
pose a cutoff limit e.g. AI = 10~^° to save integration time. 
On the other hand, if the orbit is ordered, then dj- or djj^ 
oscillate around a roughly constant mean value. This value 
is usually close to unity and in any case it is not less than 
10"^ 

Thus, by following the evolution of the smaller of the 
two indices dj- , we can distinguish between regular and 
chaotic orbits. 

We adopt a value for the Hubble time, equal to 100 half 
mass crossing times Tumct of the system. We usually calcu- 
late the LjS and AI of ordered orbits for the same number of 
radial periods, namely tj /Trj — 1200. However, in the case 
of the chaotic orbits the integration time varies due to the 
cutoff limit imposed. 

The comparison of the two methods mentioned above 
can be shown in Fig. 10. The sharp line on the left part of 
the figure (at AI^IO""*^*^) represents the most chaotic orbits 
and it is due to the cutoff limit. The triangular distribution 



on the lower right part of the figure represents the ordered 
orbits. Finally, the lane of points between these two regions 
represents orbits that are weakly chaotic. 

Almost 60% of the orbits were found chaotic. This per- 
centage is higher than the percentage of chaotic orbits in 
non rotating systems, representing elliptical galaxies, which 
is near 30% at most (see for example Voglis et al. 2002). 
Furthermore, the Lyapunov characteristic numbers of the 
chaotic orbits are larger by about one order of magnitude 
than in the non-rotating self- consistent systems, a result 
that is related to the resonance effects near corotation. How- 
ever only the 35% of the total mass can eventually develop 
chaotic diffusion in a Hubble time (Voglis et al. 2006). 

In Fig. 11 we plot the energy distribution of the parti- 
cles at the end of our calculations, that is after 250 half mass 
crossing times Thmct of the system. Figure 11a shows the dis- 
tribution of the total number of particles, where /S.Ej — 1 
(in our units which are normalized to -100 at the poten- 
tial well) and represents the amplitude of the energy bin. 
The distribution presents two peaks. Figure lib shows the 
effective potential along the X,Y and Z axes in half mass 
radius units (rhm) and Fig. 11c shows the energy distribu- 
tion of chaotic orbits alone (black line) and ordered orbits 
(gray hne). By comparing these figures we conclude that 
the sharpest peak of the distribution of the particles corre- 
sponds mostly to chaotic orbits (since there are almost no 
regular orbits there). This peak is located at Ej ^ —36, 
around corotation, as expected, and corresponds to a ra- 
dius ~ 2.2rhm. The second peak is located at Ej ^ —75, is 
related exclusively to regular orbits. 



5 STUDY OF THE ORBITAL STRUCTURE 

The purpose of this section is to classify the orbits of our 
self-consistent system in such a way as to determine the 
shape and the percentage of the orbits that support the bar. 
For this purpose we have used two powerful tools: a) the 
surface of section for test particles as well as for real N-body 
particles and b) the frequency analysis of the real N-body 
orbits. 
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5.1 Surface of section 

We fix the values of tfie coefficients of tfie analytical expan- 
sion of the potential and construct surfaces of section (SOS) 
for different energies Ej in the rotating frame (Jacobi con- 
stants) : 

Ej = ^{vl+vl + vl) + V{x, y, z) - ]^nlRly (7) 

where V(x,y,z) is the frozen potential of both the axisym- 
metric part and the bar perturbation, given by the N-body 
run, VxjVy, Vz are the velocities in the rotating frame of ref- 
erence and Qp is the angular velocity of the bar rotating as 
a solid body and calculated as in Fig. 3. 

We first use test particles for constructing SOS. In our 
3-D model the SOS are 4D and we plot their projections 
on the x,x plane, for y = 0^ y > and having initially 
z — and z — 0. All sections are made in the rotating 
frame of reference . A sim ilar technique has been u s ed by 
IContoDOulos et al.1 (|2002l ) and iKalapotharakos etHI (|200^ 1 
in order to examine the phase space structure and the fo- 
liation of the invariant tori of the orbits due to the third 
integral in non-rotating triaxial self-consistent systems. 

The values of x are normalized with the half mass ra- 
dius (hmr). The Lagrangian points Li,L2 (and therefore 
the corotation radius) are located at ~ 2.2hmr. 

In Fig. 12 we plot the SOS which corresponds to a Ja- 
cobi constant Ej = —39. All values of the Jacobi constant 
are normalized so that the well of the potential corresponds 
to -100. All the initial conditions in Fig. 12 have been cho- 
sen along x=0 and they have been integrated for 200 itera- 
tions each. Figure 12a presents the region outside corotation 
where a few stable orbits around L5 are present. Most of 
the orbits of this figure reach large distances and a part of 
them escape. Figure 12b presents the region inside corota- 
tion. Some important resonances are marked, i.e. the 3:1, 4:1 
and 5:1 resonances (that appear near the end of the bar), as 
well as the "x4" type orbit (Contopoulos and Papayanopou- 
los 1980) which has a main axis perpendicular to the bar 
and is retrograde with regard to the sign of Qp. In Fig. 13 
the corresponding periodic orbits are plotted. Every periodic 
orbit intersects the SOS at one or a finite number of points. 
The stable periodic orbits are easily identified on the SOS, 
because they are surrounded by closed invariant curves. 

In Fig. 14 we have plotted the real particles from the N- 
body run. We have chosen the initial conditions of particles 
having a Jacobi constant equal to Ej = —39 it 0.5 and we 
have integrated them in a fixed potential for 100 iterations. 
In this energy level only a small fraction of the orbits of real 
particles are regular as it can be seen from Fig. 11c. Figures 
12a and 14a present the region outside corotation containing 
test particles and real N-body orbits respectively around the 
L5 point. We see several thick dark lines that are called 
"rays". An explanation of the ray structure of this figure 
was given bv lContoDOulos Pats3 (|2Q06. ) who noticed that 
even in the chaotic region we may see a structure of rays 
(see Fig. 9 of IContopoulos k Patsis" (2006")). The successive 
iterates of every initial point along a ray are located also on 
rays, and they move in a clockwise way around the center 
of the black region. After each return the points deviate 
from their initial positions and therefore the rays thicken in 
time and finally they join each other. Thus the ray structure 



disappears if we take a much greater number of iterations 
than in Fig. 14a. 

Figure 14b presents the region inside corotation. It is 
obvious that the area corresponding to the "x4" orbits is 
depopulated while the orbits supporting the shape of the 
bar and having the same sense of rotation with it, are well 
populated. 

In Fig. 15 the surface of section for test particles is 
plotted for another value of the Jacobi constant Ej = —42. 
Figure 15a presents the region outside corotation. It is re- 
markable that some islands of stability exist that correspond 
to the OLR and -1:1 resonances. Their corresponding orbits 
are shown in Fig. 15c. For the area inside corotation (Fig. 
15b) the stable orbits correspond to the 3:1, 4:1 and 5:2 res- 
onances, as well as to the "xl" type of orbits with loops at 
the edges and the retrograde "x4". In Fig. 15d the corre- 
sponding orbits are plotted. 

The same surface of section for the real N-body par- 
ticles is shown in Fig. 16. The main resonances support- 
ing the bar are populated here again while the area around 
the "x4" orbit is depopulated. Another important remark 
is that real particles avoid the area of the OLR as well as 
the -1:1 resonance outside corotation. A remarkable feature 
common for test particles as well as for real N-body par- 
ticles is the intensely dark area around the OLR and -1:1 
resonances in the region outside corotation. (compare Figs. 
14a and 16a). This is an indication of the existence of a can- 
torus which impedes the communication of the two chaotic 
regions, imposing a partial barrier in this 2-D projection of 
the section. On the other hand since our N-body system is 3- 
dimensional, one would expect that Arnold diffusion would 
transport particles through the third dimension and pro- 
duce a communication between the two chaotic regions (for 
a definition of Arnold diffusion see Contopoulos 2002). How- 
ever, our results indicate that such a communication is very 
slow, therefore Arnold diffusion is not efficient. In Fig. 17 
the SOS is plotted for test particles, with Ej = —46 for 
the region outside corotation (Fig. 17a) and for the region 
inside corotation (Fig. 17b). Figure 17c shows the 2-D ap- 
proximation of the SOS when there is no third dimension, 
i.e. the forces on the z-axis are exactly zero. In this figure we 
mark some important tori. Every torus is characterized by 
its rotation number, i.e. the average angle between successive 
intersections of an orbit by the SOS, as seen from a central 
fixed point. This number can be represented in the form of 
a continued fraction: 

[ao,ai,...] = 5—^ — (8) 

If tti = 1 for all i beyond a certain j , this number is called a 
noble number. 

In the region outside corotation the existence of differ- 
ent tori and cantori separates the different chaotic regions. 
By comparing Fig. 17c with Fig. 17a we see only minor 
differences. In Fig. 17c the chaotic region inside the torus 
[1,1,1,7,1,1,..] is united but it is obvious that a cantorus still 
exists separating the two unequally darkened areas (a can- 
torus is a torus with a Cantor set of holes) . We have taken an 
initial condition in the most darkened area of this SOS, be- 
tween [1,1,1,7,1,1,..] and [1,2,3,1,1,..] and we have seen that 
the orbit stays located there for at least 5x10^ iterations and 
possibly for ever. Therefore we conclude that there may still 
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Figure 12. (a) The surface of section of test particles with a value of the Jacobi constant equal to -39 for the region outside corotation. 
The insert box gives the stable periodic orbit around L5 point, (b) The same for the region inside corotation. 





Figure 13. Periodic orbits that intersect the SOS at the points indicated in Fig. 12(b). 



exist tori (or cantori with very small holes) corresponding to 
a noble number [1,1,1,7,1,1..] that puts an inside barrier and 
to [1,2,3,1,1,..] that puts an outside barrier for this chaotic 
area. 

The same phenomenon is obvious in Fig. 18 where the 
same projection of SOS is plotted for real N-body particles 
for the region outside corotation (a) and for the region inside 
corotation (b). Here again it is obvious that for the region 
outside corotation, the real orbits (that are chaotic in their 
majority) are restricted in areas that are bounded by the 
same tori, while there are areas with almost no N-body par- 
ticles at all. Therefore tori and cantori on the projection of 
the SOS, play an important role in retaining chaotic orbits 
projected on the 2-D space of the galaxy, limiting the effec- 
tiveness of Arnold diffusion through the third dimension. 



The range of energies for which there exist cantori as 
described above is -60 < Ej <-45. In this range the amounts 
of ordered and chaotic orbits are of the same order (see Fig. 
11c). 

There are two mechanisms that force chaotic orbits to 
stay bound inside specific areas and make them dynami- 
cally important for the system. The first one, that has been 
described above, is the phenomenon of cantori that put a 
partial barrier on the SOS and reduce the effectiveness of 
Arnold diffusion. The second one is the stickiness effect of a 
fraction of the chaotic orbits near resonances that support 
the bar. Stickiness appears in systems of 2 or more degrees 
of freedom when chaotic orbits remain close to an invariant 
manifold for long times before escaping to lar p^e distances. 
This phenomenon is described in the book of IContopoulosI 
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Figure 14. The surface of section of real N-body particles with a value of the Jacobi constant equal to -39: (a) for the region outside 
corotation and (b) for the region inside corotation. 




Figure 15. The surface of section of test particles with a value of the Jacobi constant equal to -42: (a) for the region outside corotation 
and (b) for the region inside corotation. (c) Orbits corresponding to the region outside corotation: There are two orbits (-1:1 (b)) and a 
multiplicity 4 orbit that has bifurcated from -l:l(b). (d) Orbits corresponding to the region inside corotation. 
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Figure 16. Same as in Figure 15a,b but for the real N-body particles: (a) for the region outside corotation and (b) for the region inside 
corotation. 



(j200^). This can be shown from the frequency analysis in 
the next section. The latter mechanism makes the role of 
the chaotic layer near the limits of the bar important in 
supporting its ellipticity and boxiness. 

In Fig. 19 we plot the orbits corresponding to islands 
of stability for the region inside corotation of Fig. 17b. The 
rotation number is the average frequency of rotation around 
a fixed point of a stable periodic orbit (here OLR) on the 
phase space. The rotation number as a function of x along 
the thick line shown in Fig. 17a is plotted in Fig. 20. In this 
figure we mark several islands of stability with rational ro- 
tation numbers and the noble numbers of some important 
cantori. We n ote that the num ber of islands is smaller than 
in the case of IContopoulos ^ P atsis (2006). This is due to 
the fact that the bar perturbation is larger for the corre- 
sponding Jacobi constant in the present model, thus chaos 
is stronger. 

5.2 Frequency analysis of the orbits 

In Fig. 21 we plot the main frequencies of all the regular 
orbits out of a sample of 10764 particles at the end of the 
N-body run in polar coordinates. We have used time series 
for each orbit that correspond to a time interval of 100 ra- 
dial periods. The abscissa corresponds to the frequency of 
the variation of the radius of the orbit on the plane of rota- 
tion (x-y plane in the rotating frame of reference) and the 
ordinate corresponds to the frequency of the angle Lp be- 
tween the position vector on the plane of rotation and the 
main axis of the bar (y axis in the rotating frame of refer- 
ence). Both frequencies are regularized by the frequency in 
the third dimension (z-axis). It is obvious from this figure 
that there are groups of orbits concentrated all along Unes 
of specific resonances and others that are more scattered in 



between the lines of resonances. Moreover along the same 
Une there can be different groups of orbits belonging to the 
same resonance. For example on the 2:1 line there are three 
groups of orbits: the "xl" type with the main axis along the 
axis of the bar, the "x4" type (only a very small fraction of 
the total sample) with the main axis perpendicular to the 
main axis of the bar, and the 2:1 orbits near corotation (n.c) 
that have a more rectangular shape. 

In Fig. 22 we plot the three projections of some real N- 
body orbits that correspond to different resonances. From 
this figure it is obvious that in the majority of the orbits 
that are near resonances the third dimension is not of great 
importance and the projections on the plane of rotation have 
shapes that support the bar. We denote as "box" the orbits 
that are found scattered in between the resonances. In par- 
ticular we see some groups of orbits that are concentrated 
in rather small areas in Fig. 21, like Group A and Group 
B. Group A has orbits that support the bar. Orbits of this 
type near the limit of the bar are boxy in shape. On the 
other hand Group B orbits are more spherical on the plane 
of rotation and their projections on the other planes have a 
considerable thickness. 

Figure 23 is the same as Fig. 21 but for the chaotic 
orbits of our self-consistent system. The orbits are calcu- 
lated here up to 300 radial periods. The distribution is more 
scattered on this figure, as expected. However we can still 
identify concentrations around resonances. This is a conse- 
quence of the fact that most of these orbits are only weakly 
chaotic. 

In Fig. 24 we plot the projections of some real N-body 
chaotic orbits that have frequencies near resonances. On the 
xy-plane (of the rotating frame of reference) we have plotted 
the orbits for 100 radial periods as well as for 300 radial 
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Figure 17. SOS for Ej = —46: (a) for test particles for the region outside corotation and (b) for the region inside corotation. (c) SOS 
for the 2-D approximation of the region outside corotation. 




Figure 19. Orbits corresponding to the SOS of Fig. 17b. 

periods. We see that orbits do not change considerably from 
100 to 300 periods. 

Group A corresponds to orbits with a boxy shape and 
can be found near the end of the bar as well as in smaller 
radii. Group B corresponds to orbits more spherical and in 
general with their main axis perpendicular to the bar with 



an important thickness in the 3rd direction (z-axis). Most 
of these orbits present a "X-shape" on the zx plane. 

We have also found some families of 3-D orbits whose 
projections on the rotating plane support the bar and have 
been found and named in 3-D models by Skokos et al.l 
In Fig. 24 there are examples of real N-body chaotic 
orbits with boxy shape that support the bar all along its 
radius, like for example the 2:1 and 4:1 resonant orbits or 
the orbits of group A. Similar orbits ar e found in models 
of barred galaxies by IPatsis etiil The spatial dis- 

tribution of the chaotic orbits is in general more spherical 
than the one of the regular orbits, i.e. the shape of the bar is 
supported mostly by the regular orbits (Voglis et al. 2006). 
However there is a layer of weakly chaotic orbits that belong 
mostly to the outer part of the bar and support the shape 
of the bar. 

In Fig. 25 we show the statistics of the resonances for 
regular (gray) as well as for chaotic orbits (black) . We take a 
width ±0.02 around each resonance to populate the orbits. 
We find that 20% of the chaotic orbits lie near resonances 
and only orbits of this percentage are plotted in Fig. 24. 

In I Voglis et al.l (|2006l l the authors have found for this 



12 N. Voglis, M. Harsoula, and G. Contopoulos 




Figure 18. Same as in Fig. 17a,b for the real N-body particles. 




experiment (named QRl) that although 60% of the total 
matter is chaotic, only a fraction of 35% of the total matter 
can develop chaotic diffusion in a Hubble time. Therefore 
the rest 25% of the total matter is weakly chaotic and is 
mostly located in a layer at the edge of the bar in resonant 
form, or in regions outside corotation limited by certain tori 
and cantori (see Figs. 16,17). 



6 CONCLUSIONS 

We have studied the orbital structure of rotating self- 
consistent bar-like N-body equilibrium systems. The am- 
plitude of the bar is relatively large {Sg/go is about 80%, 
maximum). The angular velocity of rotation of the pattern 



Figure 21. Frequency analysis for the regular orbits. 

is Q p= {27r)/ {17 hmct) =0. 37 r ad /hmct, which is equivalent to 
8.2bkm/ (kpc.sec) , the spin parameter being roughly that of 
our Galaxy. 

The main conclusions of our study are the following: 

1) In our N-body experiment that simulates a rotating 
barred galaxy, chaos is very appreciable, amounting to about 
60% of the total bound matter in contrast with nonrotating 
systems, which have ^ 30% chaos. However chaos is in gen- 
eral weak and only about 35% of the mass develops chaotic 
diffusion in a Hubble time. 

2) The ellipticity of the bar decreases outwards from 0.5 
(E5 galaxy) in the inner parts to about 0.1 (El galaxy) in 
the outer parts. 
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Figure 22. Some examples of real N-body regular orbits. 
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Figure 24. Some examples of real N-body chaotic orbits. 




Figure 23. Frequency analysis for the chaotic orbits. 



Figure 25. The statistics of the regular and chaotic orbits. 



3) The bar is like an orthogonal parallelepiped with a 
boxiness parameter that is maximum in the middle of the 
bar. 

4) Only weak spirals appear in the present model, and 
they vanish relatively fast. 

5) A comparison is made between systems of test parti- 
cles and N-body systems by projecting the orbits on a 2-D 
surface of section (SOS). We notice, first, that in both cases 
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chaos is dominant outside corotation, while most orbits in 
the bar are ordered. 

6) In the outer regions of the SOS of test particles as 
well as of real N-body particles, we see some thick dark 
"rays". The particles of these rays have their images along 
other rays, but after longer times these images fill the regions 
between the rays, and the rays disappear. 

7) Certain outer regions of the SOS have different de- 
grees of darkness in the case of test particles. These regions 
are separated by invariant tori or cant or i. The invariant tori 
do not allow any communication across them in the 2-D 
case. In the 3-D case we see similar results although com- 
munication is in principle possible by means of Alnold diffu- 
sion. Thus Alnold diffusion is probably not effective in the 
present case. On the other hand cantori allow communica- 
tion across them but this process is slow if the cantori have 
small holes, and for a long time the areas from both sides of 
these cantori do not communicate. 

8) The regions inside corotation on the SOS contain 
mostly ordered orbits. However there are important differ- 
ences between the test particles case and the N-body case. In 
particular real orbits do not follow retrograde orbits around 
x4. In general N-body orbits do not enter the regions of is- 
lands if they were not inside such islands originally. 

9) We have found the forms of the periodic (planar) 
orbits and the rotation numbers along a particular straight 
line on the surface of section that allows to locate the small 
islands and the positions of tori and cantori. These tori and 
cantori have noble rotation numbers, i.e. their expansions in 
continued fractions end with an infinity of Is. 

10) A frequency analysis of ordered orbits has identified 
clearly the various resonant forms of nonperiodic orbits that 
form islands of stability. We found also stable orbits of box 
type outside the islands. 

11) We did also a frequency analysis of chaotic orbits. 
The frequency diagram has several similarities with the cor- 
responding diagram for regular orbits, although it is more 
diffuse. This is due to the fact that chaos is rather weak in 
the present case. 

12) Finally we give the statistics of the various types of 
ordered and chaotic orbits. The most important orbits are 
the xl orbits, of 2:1 type. These orbits support the bar. But 
other orbits, like the box orbits and the orbits of group A 
also support the bar. In particular the group A orbits are 
mainly responsible for the boxiness of the bar. A 20% of the 
chaotic population are located near resonances and play an 
important role in supporting the boxiness of the bar in its 
outer layers. On the other hand there is also a percentage 
of chaotic orbits that are located in restricted areas outside 
corotation limited by tori and cantori on the 2-D projections 
of SOS. The effectiveness of Arnold diffusion around these 
tori and cantori is restricted. 

In a future paper experiments with stronger and longer 
lived spiral arms will be presented. 
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